function [mimo_capacity, fig_handle] = multiuser_mimo(SNR_dB, num_users, num_antennas)
% 多用户MIMO分析
% 输入:
%   SNR_dB - 信噪比范围 (dB)
%   num_users - 用户数量数组
%   num_antennas - 天线数量数组
% 输出:
%   mimo_capacity - MIMO容量矩阵
%   fig_handle - 图形句柄

    SNR_linear = 10.^(SNR_dB/10);
    
    % 添加路径
    addpath('../Common');
    colors = color_definitions();
    
    fprintf('=== 多用户MIMO分析 ===\n');
    
    % 固定天线数量，变化用户数量
    M = 4; % 基站天线数
    K_values = num_users;
    capacity_vs_users = zeros(length(SNR_dB), length(K_values));
    
    for user_idx = 1:length(K_values)
        K = K_values(user_idx);
        fprintf('仿真 M=%d, K=%d 多用户MIMO...\n', M, K);
        
        for snr_idx = 1:length(SNR_dB)
            snr = SNR_linear(snr_idx);
            
            % 生成MIMO信道矩阵
            H = sqrt(0.5) * (randn(M, K) + 1i * randn(M, K));
            
            % 零迫波束成形
            if M >= K
                % 计算伪逆
                H_pinv = pinv(H);
                W = H_pinv' / norm(H_pinv, 'fro'); % 归一化
                
                % 计算等效信道
                H_eff = H' * W;
                
                % 计算每个用户的SINR
                for k = 1:K
                    signal_power = abs(H_eff(k,k))^2;
                    interference = sum(abs(H_eff(k,[1:k-1, k+1:K])).^2);
                    sinr = snr * signal_power / (1 + snr * interference);
                    capacity_vs_users(snr_idx, user_idx) = capacity_vs_users(snr_idx, user_idx) + log2(1 + sinr);
                end
            else
                % 天线数少于用户数，使用时分多址
                for k = 1:K
                    user_channel = H(:,k);
                    channel_gain = norm(user_channel)^2;
                    user_capacity = (K/M) * log2(1 + snr * channel_gain / K);
                    capacity_vs_users(snr_idx, user_idx) = capacity_vs_users(snr_idx, user_idx) + user_capacity;
                end
            end
        end
    end
    
    % 固定用户数量，变化天线数量
    K = 4; % 用户数
    M_values = num_antennas;
    capacity_vs_antennas = zeros(length(SNR_dB), length(M_values));
    
    for ant_idx = 1:length(M_values)
        M = M_values(ant_idx);
        fprintf('仿真 K=%d, M=%d 多用户MIMO...\n', K, M);
        
        for snr_idx = 1:length(SNR_dB)
            snr = SNR_linear(snr_idx);
            
            % 生成MIMO信道矩阵
            H = sqrt(0.5) * (randn(M, K) + 1i * randn(M, K));
            
            % 使用不同预编码方案
            if M >= K
                % 零迫波束成形
                H_pinv = pinv(H);
                W_zf = H_pinv' / norm(H_pinv, 'fro');
                
                % MMSE波束成形
                W_mmse = H * inv(H' * H + (K/snr) * eye(K)) / sqrt(K);
                
                % 计算和速率
                for scheme = 1:2
                    if scheme == 1
                        W = W_zf;
                    else
                        W = W_mmse;
                    end
                    
                    H_eff = H' * W;
                    sum_rate = 0;
                    
                    for k = 1:K
                        signal_power = abs(H_eff(k,k))^2;
                        interference = sum(abs(H_eff(k,[1:k-1, k+1:K])).^2);
                        if scheme == 1
                            sinr = snr * signal_power / (1 + snr * interference);
                        else
                            noise_power = 1 + snr * interference;
                            sinr = snr * signal_power / noise_power;
                        end
                        sum_rate = sum_rate + log2(1 + sinr);
                    end
                    
                    if scheme == 1
                        capacity_vs_antennas(snr_idx, ant_idx) = sum_rate;
                    end
                end
            else
                % 天线数不足，使用时分多址
                sum_rate = 0;
                for k = 1:K
                    user_channel = H(:,k);
                    channel_gain = norm(user_channel)^2;
                    user_capacity = (K/M) * log2(1 + snr * channel_gain / K);
                    sum_rate = sum_rate + user_capacity;
                end
                capacity_vs_antennas(snr_idx, ant_idx) = sum_rate;
            end
        end
    end
    
    % 绘制结果
    fig_handle = figure('Name', '多用户MIMO分析', 'Position', [250, 250, 1200, 800]);
    
    % 容量 vs 用户数量
    subplot(2,4,1);
    for user_idx = 1:length(K_values)
        plot(SNR_dB, capacity_vs_users(:, user_idx), ...
             ['-', colors(user_idx)], 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('SNR (dB)');
    ylabel('和速率 (bps/Hz)');
    title('容量 vs 用户数量 (M=4)');
    legend(arrayfun(@(x) sprintf('K=%d', x), K_values, 'UniformOutput', false));
    
    % 容量 vs 天线数量
    subplot(2,4,2);
    for ant_idx = 1:length(M_values)
        plot(SNR_dB, capacity_vs_antennas(:, ant_idx), ...
             ['-', colors(ant_idx)], 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('SNR (dB)');
    ylabel('和速率 (bps/Hz)');
    title('容量 vs 天线数量 (K=4)');
    legend(arrayfun(@(x) sprintf('M=%d', x), M_values, 'UniformOutput', false));
    
    % 用户数量对频谱效率的影响
    subplot(2,4,3);
    spectral_efficiency = capacity_vs_users ./ (K_values' * ones(1, length(SNR_dB)))';
    for user_idx = 1:length(K_values)
        plot(SNR_dB, spectral_efficiency(:, user_idx), ...
             ['-', colors(user_idx)], 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('SNR (dB)');
    ylabel('频谱效率 (bps/Hz/用户)');
    title('每用户频谱效率');
    
    % 天线效率
    subplot(2,4,4);
    antenna_efficiency = capacity_vs_antennas ./ (M_values' * ones(1, length(SNR_dB)))';
    for ant_idx = 1:length(M_values)
        plot(SNR_dB, antenna_efficiency(:, ant_idx), ...
             ['-', colors(ant_idx)], 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('SNR (dB)');
    ylabel('天线效率 (bps/Hz/天线)');
    title('每天线效率');
    
    % 不同预编码方案比较
    subplot(2,4,5);
    M = 4; K = 4; snr = 15; % dB
    snr_linear = 10^(snr/10);
    
    % 生成信道
    H = sqrt(0.5) * (randn(M, K) + 1i * randn(M, K));
    
    % 零迫波束成形
    H_pinv = pinv(H);
    W_zf = H_pinv' / norm(H_pinv, 'fro');
    H_eff_zf = H' * W_zf;
    
    % MMSE波束成形
    W_mmse = H * inv(H' * H + (K/snr_linear) * eye(K)) / sqrt(K);
    H_eff_mmse = H' * W_mmse;
    
    % 计算用户速率
    user_rates_zf = zeros(1, K);
    user_rates_mmse = zeros(1, K);
    
    for k = 1:K
        % ZF
        signal_power = abs(H_eff_zf(k,k))^2;
        interference = sum(abs(H_eff_zf(k,[1:k-1, k+1:K])).^2);
        sinr_zf = snr_linear * signal_power / (1 + snr_linear * interference);
        user_rates_zf(k) = log2(1 + sinr_zf);
        
        % MMSE
        signal_power = abs(H_eff_mmse(k,k))^2;
        interference = sum(abs(H_eff_mmse(k,[1:k-1, k+1:K])).^2);
        sinr_mmse = snr_linear * signal_power / (1 + snr_linear * interference);
        user_rates_mmse(k) = log2(1 + sinr_mmse);
    end
    
    bar(1:K, [user_rates_zf; user_rates_mmse]');
    grid on;
    xlabel('用户索引');
    ylabel('速率 (bps/Hz)');
    title('预编码方案比较 (SNR=15 dB)');
    legend('ZF', 'MMSE');
    set(gca, 'XTick', 1:K);
    
    % 信道条件数分析
    subplot(2,4,6);
    condition_numbers = zeros(1, 1000);
    for i = 1:1000
        H = sqrt(0.5) * (randn(M, K) + 1i * randn(M, K));
        condition_numbers(i) = cond(H);
    end
    histogram(condition_numbers, 50, 'Normalization', 'pdf');
    grid on;
    xlabel('条件数');
    ylabel('概率密度');
    title(sprintf('信道矩阵条件数分布 (M=%d, K=%d)', M, K));
    
    % 容量缩放律
    subplot(2,4,7);
    K_range = 2:2:16;
    M_fixed = 8;
    scaling_rates = zeros(1, length(K_range));
    snr = 15; % dB
    snr_linear = 10^(snr/10);
    
    for k_idx = 1:length(K_range)
        K = K_range(k_idx);
        H = sqrt(0.5) * (randn(M_fixed, K) + 1i * randn(M_fixed, K));
        
        if M_fixed >= K
            H_pinv = pinv(H);
            W = H_pinv' / norm(H_pinv, 'fro');
            H_eff = H' * W;
            
            sum_rate = 0;
            for k = 1:K
                signal_power = abs(H_eff(k,k))^2;
                interference = sum(abs(H_eff(k,[1:k-1, k+1:K])).^2);
                sinr = snr_linear * signal_power / (1 + snr_linear * interference);
                sum_rate = sum_rate + log2(1 + sinr);
            end
            scaling_rates(k_idx) = sum_rate;
        else
            sum_rate = 0;
            for k = 1:K
                user_channel = H(:,k);
                channel_gain = norm(user_channel)^2;
                user_capacity = (K/M_fixed) * log2(1 + snr_linear * channel_gain / K);
                sum_rate = sum_rate + user_capacity;
            end
            scaling_rates(k_idx) = sum_rate;
        end
    end
    
    plot(K_range, scaling_rates, 'o-', 'LineWidth', 2);
    hold on;
    plot(K_range, min(K_range, M_fixed) * log2(1 + snr_linear), 'r--', 'LineWidth', 2);
    grid on;
    xlabel('用户数量');
    ylabel('和速率 (bps/Hz)');
    title('容量缩放律');
    legend('仿真结果', '理论上限');
    
    % 能量效率
    subplot(2,4,8);
    power_range = 5:5:30; % dBm
    energy_efficiency = zeros(1, length(power_range));
    K = 4; M = 4;
    
    for p_idx = 1:length(power_range)
        power = 10^((power_range(p_idx) - 30)/10); % 转换为瓦特
        snr_linear = power / 1e-3; % 假设噪声功率为-30dBm
        
        H = sqrt(0.5) * (randn(M, K) + 1i * randn(M, K));
        H_pinv = pinv(H);
        W = H_pinv' / norm(H_pinv, 'fro');
        H_eff = H' * W;
        
        sum_rate = 0;
        for k = 1:K
            signal_power = abs(H_eff(k,k))^2;
            interference = sum(abs(H_eff(k,[1:k-1, k+1:K])).^2);
            sinr = snr_linear * signal_power / (1 + snr_linear * interference);
            sum_rate = sum_rate + log2(1 + sinr);
        end
        
        energy_efficiency(p_idx) = sum_rate / power; % bits/Joule
    end
    
    plot(power_range, energy_efficiency, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('发射功率 (dBm)');
    ylabel('能量效率 (bits/Joule)');
    title('能量效率分析');
end